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ABSTRACT 

We present a general expression for a lognormal filter given an arbitrary nonlinear galaxy 
bias. We derive this filter as the maximum a posteriori solution assuming a lognormal prior 
distribution for the matter field with a given mean field and modeling the observed galaxy dis- 
tribution by a Poissonian process. We have performed a three-dimensional implementation 
of this filter with a very efficient Newton-Krylov inversion scheme. Furthermore, we have 
tested it with a dark matter N-body simulation assuming a unit galaxy bias relation and com- 
pared the results with previous density field estimators like the inverse weighting scheme and 
Wiener filtering. Our results show good agreement with the underlying dark matter field for 
overdensities even above 5 ~ 1000 which exceeds by one order of magnitude the regime in 
which the lognormal is expected to be valid. The reason is that for our filter the lognormal 
assumption enters as a prior distribution function, but the maximum a posteriori solution is 
also conditioned on the data. We find that the lognormal filter is superior to the previous filter- 
ing schemes in terms of higher correlation coefficients and smaller Euclidean distances to the 
underlying matter field. We also show how it is able to recover the positive tail of the matter 
density field distribution for a unit bias relation down to scales of about > 2 Mpc/h. 

Key words: (cosmology:) large-scale structure of Universe - galaxies: clusters: general - 
catalogues - galaxies: statistics 



1 INTRODUCTION 

The luminous matter we observe on the sky represents only a small 
fraction of the total matter in the Universe and yet with a careful 
treatment of the observational selection effects and the processes 
of galaxy formation we can hope to extract valuable information 
about the distribution of all matter from the distribution of lumi- 
nous matter alone. The more precise the techniques for making this 
connection are the better we will be able to test our theories for the 
history of the Universe. 

In 1934 Hubble found that the distribution of galaxy c ounts in 
cells on the sky is we ll fitted by a logn ormal distribution dHubbld 
Il934l) . More recently. Iwild et al.1 12005) showed that this model is 
valid at least down to gridding scales of about 10 Mpc for galax- 
ies in the 2DF catalogue. As galaxies are good tracers of matter 
on large cosmological scales the lognormal model should also ap- 
ply to the matter field at least to some degree. Kita ura etaD J2009) 
showed recently that the matter field reconstruction using (least 
squares) Wiener filtering is very well fit by a lognormal distribu- 
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tion after smoothing the reconstruction with a Gaussian kernel of 
radius rs for 10 Mpc < rs < 30 Mpc. 

From a physical point of view, one would expect the density 
field to be lognormally distributed after it has been smoothed on an 
appropriate scale. This follows from assuming an initially Gaussian 
density and velocity field and extrapolating the continuity equation 
for the matter flow into the nonlinear reg ime with the linear veloc- 
ity fluctuations (see lColes & Jonesll99ll) . Since the lognormal field 
is not able to describe caustics, we exp ect this dis t ribut ion to fail 
below some threshold smoothing scale. iKavo et all d200lh demon- 
strated that the lognormal distribution is a good approximation up 
to overdensties of about S ^100. 

Shortly after the succes s of Wiener filterin g in large-scale 
structure reconstruction (see IZaroubi et al which assumes 

a Gaussian prior for the matter field, a reconstructio n filter base d 
on the lognormal prior distribution was proposed fsee lShethll 19951) . 
With such a filter nonlinea rities in the density field should be better 
recovered. In[sheth (1995) the Wiener filter was generalized to be 
applied to a lognormal distribution by a variable transformation. A 
problem with this approach is that the noise covariance has a com- 
plex form even for the simple Poisson likelihood assumption and is 
difficult to efficiently apply to realistic data-sets. 
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The idea of modeling the galaxies as Poisson-sampled 
from a lognormal un derlying field was first applied to data by 
Saunders &Ballinger (2000). They proposed to use a filtering 
scheme based on an expansion of the logarithm of the matter field 
as a sum of harmonics . The density reconstru ction using this tech- 
nique as presented by Saunders et al. (2000) is fairly smooth and 
nonlinear structures cannot be easily recognised. This could be due 
to the sparse sampling of the PSCz catalogue which was used in 
their study or to the trunca tion of the harm onic series. 

As demonstrated in lKitaura~& EnBlin (2008), the Poissonian 
likelihood can be easily regularised by combining it with a prior 
when estimating the maximum a posteriori. They showed this cal- 
culation for Gaussian and entropic prior distribution functions. 

The idea of using the full Poissonian likelihood without re- 
maining at second order approximations using only the noise co- 
variance mat rix is based on the Richardson-Lucy deconvolu t ion al - 
gorithm fsee lRichards^ll97llLu^ll974h . [shepp & Vardil dl982h 
showed that this filter comes from the maximum likelihood esti- 

1 J ~T 

Nusser & Haehnelt ( 1999) pro- 
posed using this method to recover the density field from the Ly- 
man alpha forest. The problem that arose here was that the algo- 
rithm requires truncation as it assumes a flat prior for the matter 
field and thus the decon volution of the respons e operator is not 
regularised. However, as Kitaura & EnBlin (2008) point ed out, this 
kind o f problem can be solved by introducing a prior. lEnBlin et al.l 
(2008) proposed to calculate higher order corrections to obtain an 
estimate for the mean of the posterior distribution by employing a 
generating functional formalism with a Poissonian process on top 
of a lognormal field for the galaxy distribution. 

In this work we present a general expression for the Poisson- 
lognormal filter given an arbitrary nonlinear galaxy bias. We derive 
this filter as the maximum a posteriori solution assuming a log- 
normal prior distribution for the matter field with a constant mean 
field and modeling the observed galaxy distribution by a Poisso- 
nian process. We have performed a three-dimensional implemen- 
tation of this filter with a very efficient Newton-Krylov inversion 
scheme extending the ARGO computer code to perform nonlin- 
ear inversions (see IKitaura & E nBlin 2008). Furthermore, we have 
tested it for a linear galaxy bias relation and compared the results 
with other density field estimators commonly used in the litera- 
ture (e.g. the inverse weighting scheme and the least squares (LSQ) 
Wiener filter). The one-dimensional lognormal probability distri- 
bution is known to fit the matter distrib ution wel l up to overden- 
sities of about 5 ^100 as found by iKavo et all (1200 ll) . Our re- 
sults show, however, good agreement for overdensities even above 
S ^1000 which exceeds by one order of magnitude the expected 
regime in which the lognormal is expected to be valid. The reason 
for this apparent disagreement is that for the filter presented here 
the lognormal assumption enters as a prior distribution function, 
but the maximum a poster iori solution i s also conditioned on the 
data. For the same reason iKitaura et al.l J2009h obtained a highly 
non-Gaussian distributed matter field after using LSQ-Wiener fil- 
tering which according to the Bayesian formalism assumes a Gaus- 
sian distribution. We find that the Poisson-lognormal filter has a 
range of applicability in recovering matter density fields down to 
scales of about > 2 Mpc/h. However, the matter statistics show that 
the Poisson-lognormal filter fails to recover underdense regions 
5 < — 0.6 with very few data. In addition, we test the maximum a 
posteriori assuming a Gaussian prior and found that it is not capa- 
ble of recovering the density field when 5 ^> 1 and gives negative 
densities in low density regions which makes this filter unreliable 
for recovering densities of 5 < 1. 
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Figure 1. Two models of completeness w emulating apparent magnitude 
limit effects (continuous curve: w\ and dashed curve: W2) dependent on 
the distance r to the observer in Mpc/h. 

Finally, we show in appendix [A] that the LSQ filter is the op- 
timal linear filter under a Poisson noise assumption using up to 
second order statistics and does not neglect any signal to noise cor- 
relation, co ntrary to what has been assumed in the literature (see 
for examp le Zaroub i_et al l l 1 9951 : ISeli ald[l998l : Erdog du et al . 2004; 
IKitaura & EnBlinl2008l) . We also derive in appendix|B]a filter with a 
lognormal model for the underlying signal and an additive, signal- 
independent and Gaussian distributed noise which could be of in- 
terested in other fields of astronomy. 

The paper is structured as follows. In section we present the 
Bayesian approach used in this work. After defining the likelihood 
for the galaxy sample and the prior distributions for the matter field 
we calculate the maximum a posteriori (MAP) estimates for the un- 
derlying density field. Then the numerical scheme is presented in 
section [3] which permits us to solve the MAP estimates. We then 
present in section |4] a series of numerical experiments which show 
the performance of the different density estimators. Finally, we dis- 
cuss our results. 



2 BAYESIAN APPROACH 

A Bayesian approach requires the definition of a likelihood and a 
prior. A full Bayesian analysis would require the complete char- 
acterizati on of the posterior d istribution using sampling schemes 
(see e.g. IWandelt et al .1 12004). We leave such an approach for a 
forthcoming publication and restrict ourselves here to calculate the 
extrema which leads to the maximum a posteriori expressions. This 
permits us to get a fast estimate of the density field. In this work, 
we consider a Poissonian likelihood for the observed distribution 
of galaxies and combine it with a Gaussian and a lognormal prior 
distribution for the overdensity field. In the next subsections these 
distribution functions are presented and the calculation of the dif- 
ferent MAP-estimators are shown in detail. 

2.1 Poissonian likelihood 

The likelihood represents the observation process which leads to 
the data. It is the probability distribution function that describes 
the nature of the observable. In this case we look for a model that 
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accounts for the discrete nature of a galaxy distribution, the so- 
called shot noise. This kind of noise is t raditionally m odeled by 
a Poissonian distribution (see for example Peebles 1980a). Such a 
model assumes that each cell of the Universe in which we count 
some number of galaxies (maybe according to a certain luminos- 
ity type) is statistically independent from each other. However, the 
variance of counts in cells including a correlation term predicts 
the non-Poisso nian character of the distribution of galaxies (see 
Peebles 1980a). Hierarchical structure formation models assume 
that galaxies form inside d ark matter halos via t h e energy dissipa- 
tion b y baryons (see e. g. I White & Reeslll978l) . ISomerville et al.l 
booih showed based on numerical N-body simulations that in re- 
gions of lower than average overdensity, the scatter in the halo bias- 
ing (the relation between the dark matter halos and the underlying 
dark matter distribution) is generally smaller than the mean Poisson 
shot noise, and in overdense regions i t is larger (for a revie w on the 
halo model see lCoorav & Shethl2002h . lMo & White! dl996h already 
pointed out that halo-exclu sion can cause sub-Poisson variance. 
ICasas-Miranda et all d2002h demonstrated with higher resolved N- 
body simulations that the galaxy biasing process, as well as the 
halo biasing process, is not only determined by the local value of 
the mass density field, but also by other local quantities, such as 
dumpiness, and by non-local properties, such as the large-scale 
tidal field. Accounting for all these effects is out of scope of this 
work, but should certainly be further investigated. 

Here, we will restrict ourselves to a model in which the ob- 
served distribution of galaxies is given by an inhomogeneous Pois- 
son realization of a c ontinuous den sity field. We define the likeli- 
hood function as (seeEeebles 198q3): 



C((N° g ) g \N° g ) = Htt ls exp [-(N^), 



<^>g 



N2A 



(1) 



with N gji denoting the number count of observed galaxies in cell 
i, and N ce \\ s being the total number of cells. Here ({ }) g = 
({ })(jv||a°) = Div°=o Ppois(N° | w\){ } denotes an ensemble 
average over the Poissonian distribution with the expected num- 
ber of galaxy counts given by the Poissonian ensemble average: 
A° = wX = (Ng) s . The expected number count is related to the 
underlying continuous galaxy overdensity field S gl i through: 



<AT g ° ;i ) g = iV g ^(l + J g ,0, 



(2) 



where N g is the mean number count of galaxies and w% the com- 
pleteness at cell i. The logarithm of the likelihood can be written 
as: 

ln£ 2 = -iV g ^(l + ^ g ,0+K^ n (^g^(l + ^))-ln(^ g °z!). 

(3) 



2.2 Gaussian prior 

The prior probability distribution function describes the statistical 
nature of the signal one wants to infer from the observed data. 
Here the physical model of the underlying matter field comes 
in. As inflationary scenarios predict a close to Gaussia n distribu- 
tion function for the initial density fluctuations (see iGuthl 1 1 98 1 |: 
Guth & Pilll982l: IStarobinskyl Il982l: hawking! 1 19811 : lLinddll982l : 
Albre cht & Steinhardtlll982l : lBardeen et al.ll983l) "and linear theory 
preserves this property throughout cosmic evolution it is reasonable 
to assume a Gaussian prior to model the large-scale matter field. 
Note, however, that this can only be true for |<5| <C 1 since oth- 
erwise the Gaussian di stribution predicts unp hysical negative den- 
sities. Here we follow Bardee net al.l dl986h to describe the prior 



probability distribution of the density field by a multivariate Gaus- 
sian distribution function: 



V(S M \p) 



V^TT^ceiis det(S) 



exp 



(4) 



with p being the set of cosmological parameters which determine 
the autocorrelation matrix S and 5m is the overdensity in mass. 
The application of the autocorrelation matrix S to a vector x is a 
convolution of the form: Sx = S(r) o x(r) (with r being the cell 
coordinates in configuration space of the box r = {ix ,iy,iz} and 
with "o" denoting the convolution operation). 

Note that the Fourier transform of the autocorrelation matrix S 
is equal to the power spectrum: S(fe, k') = (27r) 3 P(fe / )fo(fe — fc' ) 
(using the same Fourier definitions as in iKitaura & EnB lin 2008). 
The logarithm of the prior distribution function can be written as: 

\nV(5 M \p) = -i^S-^M + c, (5) 

with c being the logarithm of the normalization. 

The posterior distribution function P is proportional to the 
product of the prior V and the likelihood C. To find the maximum 
a posteriori (MAP) we need to calculate the extremum. Performing 
the derivative of the posterior with respect to the matter overdensity 
field 5m yields: 



0\nP 0\nV 0\nC n 
oc —777; h — — 0. 



05 M 05 M 
The derivative of the prior leads to: 
d\nV 



05* 



05* 



(6) 



(V) 



Since the likelihood is expressed as a function of the galaxy density 
field, we need to define the bias between the galaxy and matter 
fields. 



2.2.7 Linear bias 

As a particular case, let us consider a linear bias function given by: 

5g,i = bjj5M,ji (8) 
3 

which relates the corresponding power spectra in the following 
way: b(k) — ^/P g (k) / Pwi(k), with P g (k) being the galaxy 
power-spectrum and Pyi(k) being the matter power-spectrum. 

The derivative of the likelihood with respect to the matter 
overdensity fields is then given by: 

N° 

-N g Wi + 



1 + J2i bi,iS M ,i_ 



(9) 



Adding this result to the prior term Eq. we obtain the MAP 
equation: 



3 1 v 



N° 1 



N g wi 



(10) 



with the superscript G standing for the Gaussian prior assumption. 



2.2.2 Unity bias 

Let us consider the special case when the matter field is equal to a 
continuous galaxy field: 



5x.i — 5^ 



(id 



© 0000 RAS, MNRAS 000, 000-000 



4 Kitaura et al. 



then the MAP equation reads: 



(12) 



2.3 Lognormal prior 

Here we introduce th e lognormal prior distribution as proposed by 
IColes& JonesUl99ll) : 



V{s\p) 



y/(2ir) N ceiis det(S L 



: exp 



(13) 



with s being the logarithm of the weighted matter density: 

s ^ = log(l + <5 M i), (14) 

and Sl the corresponding autocorrelation matrix. Note that the log- 
normal autocorrelation function Sl applied to a vector x is again a 
convolution: Slx = Sh(r)ox(r). The transformation of the corre- 
lation function corresponding to the overdensity field to the signal 
s is given by: 

5 L (r) = Iog(l + 5(r)). (15) 

The mean field fi is taken to be: 

= -erg/2, (16) 

with erg = Sl(0) as used by iKavo et al.1 d2QQlh to ensure an 
overdensity field w ith zero mearQ (for a formal derivation see 
IColes & Jo nes 1991). The logarithm of the prior distribution yields: 

\nV(s\p) =-\{s- tfSZ\a - m) + c, (17) 

with c being some constant term. In this case, we look for the ex- 
tremum with respect to the signal s: 



d\nP 



d\nV d\nC 

+ 



0. 



(18) 



ds ds ds 

The derivative of the prior has now an additional term due to the 
mean-field /x: 

dlnV 



ds 



-Sl 1 (a - /x) • 



(19) 



Now we need to relate the galaxy field to the matter field in or- 
der to express the likelihood as a function of the signal statistically 
defined through the prior distribution function. 



The factor relating the galaxy field to the matter field yields: 

ain(l + * s ,0 = ain(l + B(exp(*)-r) t ) 
d\n(l + 5 M ,k) ds k 

dB(8 M )i (1 + M 



(23) 



dS M ,k 1 + B(6m)i' 

The final result for the derivative of the likelihood with respect to 
the logarithm of the normalized matter field s assuming a general 
nonlinear bias is given by: 

d\nd 



E : 



ds k 



(24) 



EdB(S M ) i (l+#M,Jfe) / ITT /-, . d /jj \ \ I AT o \ 
i i tmx \ {-NgWi 1 + B(S M )i) + N Sii ) . 
oo M ,k 1 + B(d M )i 

i 

Combining this result with the prior term (Eq. \\9l we obtain the 
MAP equation: 

^S^(ln(l+4,,)-w) = 



(25) 



dB(6l l ) l (1+6^ 
, dS^ 1 + 3(61, 



E 



with the superscript L standing for the lognormal prior. 



2.4.1 Linear bias 

For the linear bias case the derivative of the likelihood reduces to 
the following expression: 



E 



<9 In Cj _ ^ &i,fc(l + 6 M ,k) 
ds k A" 1 + 52 v b itl t6 Mtl t 



(26) 



Accordingly, the MAP equation reads: 
^5^(ln(l + ^)- W ) = 



(27) 



2.4 General nonlinear bias 

Let us consider here a general nonlinear relation between the galaxy 
field and the matter field. 



B(6* 



(20) 



The derivative of the likelihood with respect to the signal s which 
we want to recover can be written as: 



Ed In Hi 



din d dln(l + a g ,Q 
^ ain(l + (J g ,z) ain(l + (5 M ,fc)" 



(21) 



The derivative of the likelihood with respect to the logarithm of the 
normalized galaxy field ln(l + 6 gi i) yields: 

<91n d 



<91n(l + £ g;Z 



(iY g ° ; , - N s wi (1 + B(5 M )i)) 6%. (22) 



1 Here generalized to a multivariate lognormal distribution. 



2.4.2 Unity bias 

For an unity bias the derivative of the likelihood reduces to: 

S g ,i = 6m, i- (28) 
Then, the derivative of the likelihood reduces to: 



E 



dlnCi 



-N g w k exp (s k ) + iVg jfc , (29) 



ds k 

and the MAP equation reads: 

^2 S ^lo ( s i ~ Vj) = N li ~ N s w i exp(si). (30) 

3 

Using the definitions: N g wj exp(sj) = N s Wj(l + 6 gj i) — (N°) g 
and e°j = N° — (N°) g we can rewrite the MAP equation as: 



(31) 
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The signal s is thus given by the propagation of the noise e° given 
the correlation Sl up to a shift due to the mean fi. Expressing it as 
a function of the matter overdensity-field we obtain: 

E ( ln {} + *M,i) - to) = (Ki ~ N e wi(l + . 

(32) 

3 NUMERICAL APPROACH 

The problem we are studying here requires the solution of a nonlin- 
ear system of 256 3 (about 17 • 10 6 ) coupled equations. Note, that 
each cell introduces an equation. Thus, to find the MAP solution 
(Eqns. Q2] and [26]) we apply an operator base d iterative inversion 
scheme as proposed in lKitaura & EnBlinl (120081) which reduces the 
most expensive operations to FFTs. In particular, we use a nonlin- 
ear Newton-Krylov scheme which i s briefly presented in th e next 
subsections (for a reference see e.g. iKitaura & EnBlinl l2008l) . 

3.1 Method 

Let us write a system of nonlinear equations as: A(x) — /, with A 
being the nonlinear operator dependent on x and / some constant 
vector. We then define the gradient of the quadratic approximation 
as: 

VQ(x) = A(x) - f. (33) 

The corresponding Hessian matrix is then given by the second 
derivative of the gradient of Q: 

H = VVQ(x). (34) 

The basic Newton-Raphson solver scheme is given by: 



hl x j - (H j ) 1 VQ{x j 



(35) 



This scheme turns out to be extremely inefficient. Therefore, we 
implement a Krylov step in which the solution is updated in the 
following way: 

x^ 1 = X J +r j e, (36) 



with the stepsize r given by (for a derivation see lKitaura & EnBlinl 
l2008h : 



^vg(a! J 



(37) 



and £ being the search ing vector (for schemes to calculate £ see 
IKitaura & EnBlin 2008). In the next subsections we give the par- 
ticular expressions for the quantities required to calculate the MAP 
given a Gaussian prior first and finally given a lognormal prior. 



3.1.1 Gaussian prior 

Rewriting Eqn.[T2las: 

<*S -Sdiag(l + (5^)- 1 iV° = N g Sw g . (38) 

We can identify A(s) = - Sdiag(l + (5^) _1 iV° and / = 
NgSw g . The corresponding gradient of the quadratic form is given 
by: 



VQ(s) = <5m - S (diag(l + ^m)" 1 ^ - N s w f 
and the Hessian yields: 

H = l + Sdiag(l + ^)- 2 AT°. 



) , (39) 



(40) 



3.1.2 Lognormal prior 

We formulate Eq.[32]in an analogous way to the previous subsec- 
tion as: 

1 (* - M) + diag (~N g w g ) exp(s) = N° s , (41) 

with A(s) = S~ 1 (s - ii) + diag (~N g w g ) exp(s) and / = JV°. 
Thus, the gradient of the quadratic form is given by: 

VQ(s) = Sl 1 (a - /i) + diag (iV g ^ g ) exp(s) - N° e , (42) 

and the corresponding Hessian matrix reads: 

H = Sl 1 + diag (N g w g ) diag (exp(s)) . (43) 
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LIKELIHOODS 

Gaussian Poissonian 



PRIORS 

Flat (a) 5™ = -5k - 1 

w^Ng 

Gaussian = £ . (s~j + wjN^)' 1 7V g 5° . . = E, Si,,- ( J3g_ - JV gWj 

Lognormal (b) £\ S^. (in (l + ifc J _ ^) = jv° . _ N gWi (l + ^ .) 

Table 1. Filters which are used in this work classified by the assumed likelihood and prior (with the exception of (a) and (b)). Note, that the bias has been set 
to one. (a) COBE-filter used in CMB mapping (see Janssen & Gulkis 1992). (b) for a derivation of this filter see appendix [B] 



4 NUMERICAL EXPERIMENTS 

In this section we investigate the performance of the Poisson- 
lognormal filter. We construct the mock observed galaxy distribu- 
tions by making a Poisson sample over the dark matter particles 
of the Millennium run according to different completeness models 
(see subsection 14.21 ). This permits us to avoid the galaxy biasing 
and redshift distortions problems in our tests. 

We also test the Poisson-lognormal filter against other filters. 
We make a comparison to the MAP with a Gaussian prior assump- 
tion, to the inverse weighted data, and to the LSQ-Wiener filter 
(see section below l47TTTb . For an overview of the filters used in this 
work see table [T] 

Finally, we test the quality of the reconstruction by making 
a cell-to-cell comparison to the underlying matter field which is 
assumed to be given by the da rk matter distributio n of the Millen- 
nium run at redshift zero fsee lSpringelet al. 2005). In addition, we 
study the matter density statistics of the dark matter field and the 
reconstructions. 



4.1 Quality validation of the density reconstruction 

In order to show the performance of the Poisson-lognormal fil- 
ter we compare the results with two othe r est imates of the density 
field. We follow lKitaura&EnBlinl (120081) and lKitaura etaD <2009t) 
in quantitatively measuring the quality of the reconstructions. 



4.1.1 Alternative density field estimators: Inverse weighting and 
LSQ-Wiener filtering 

Let us first introduce a representation of the data which tries to 
compensate for the selection function effect which we call inverse 
weighting (IW). We define the inverse weighted galaxy number 
count per cell i as: 



iV IW = — N°- 



(44) 



The corresponding inverse weighted overdensity is calculated as 
follows: 



N 



elVV _ _ 'g,2 



(45) 



Note that the inverse weighting scheme can be derived as the max- 
imum likelihoo d estimator assumin g a Poissonian li kelihood (for 
a derivation see Kita ura et alj l2009). As discussed in Kitaur a et al. 1 
(2009) IW boosts the estimated density field at low completeness. 
Therefore it includes in gen eral an additional sm oothing step which 
lessons this effect (see e.g. lErdogdu et alj|2004l) . 

For an additional comparison let us introduce the least squares 



version of th e Wiener filter (or LSQ filter for short) given by 
iKitaura et al.1 (see l2009l) : 

(5m SQ = B" 1 (s _1 + W+N^w) _1 W^" 1 ^, (46) 

with W being the three dimensional mask operator defined by: 
Wij = WjS^j (Sf^j is the Kroenecker delta), and the Fourier trans- 
form of B given by: B kjk > = by5^ k , as introduced in subsection 
12.4.11 We define the observed galaxy overdensity which we use as 
the input vector for the LSQ reconstruction by: 



N° 

g ' 2 - N s 



(47) 



The n oise term in Eq.[46]has the following form fsee lKitaura et al] 
Hqq3): 



Wi K 



(48) 



Note that the LSQ-Wiener filter is the optimal linear filter using 
up to second order statistics even for the Poisson-noise assumption 
for which the noise is signal dependent. There is not an additional 
assumption or approximation by neglecting the signal to noise cor- 
relatio n. This point has been unclear in the literature (see fo r ex- 
ample [z^oubLet^li^; [1^ We 
show that the signal and the noise are indeed uncorrected in the 
appendix. The LSQ-Wiener filter also happens to be the MAP fil- 
ter for a Gaussian likelihood and a Gaussian prior as indicated in 
tabled] 

In our numerical experiments we use a unity bias: b k = 1. 
Thus, the Fourier transform of S yields: P g (k) = Pm(Aj). The 
power-spectrum Pyi(k) is given by a nonlinear fit that also de- 
scrib es the effects of virialised structures with a halo term as given 
bv ISmithetalJ (l2003h at redshift z = 0. We choose the con- 
cordance A CDM-cosmology with Q m — 0.24, Qk — and 
= 0.76 dSpergel et alj2007b . In addition, we assumed a Hubble 
constant with h = 73 and a spectral index n s = 1. 



4.1.2 Quantitative measures 

Let us define the correlation coefficient r between the reconstructed 
and original matter density fields by: 

Ef cells Mr 



r(<r c ,j M ) 



V£f cells (<5M,i)V£f cellB (* 

and the Euclidean distance 



(49) 



rec \ 
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4.2 Input data setup 

We construct the mock observed galaxy distribution taking a ran- 
dom subs ample of the particl es in the Millennium run at redshift 
zero fsee lSpringel et all2005h which was gridded on a 256 3 mesh. 
Later we also investigate the resolution dependence using a 128 3 , 
and a 64 3 mesh. As already stated above, our setup permits us to 
avoid the galaxy biasing problem in our tests. Note, that we also 
avoid the redshift distortions by considering the dark matter par- 
ticles in real-space. In this way we generate three different input 
mock galaxy catalogues. One has about 1 Million particles and has 
been produced as a Poisson sampling with a homogeneous com- 
pleteness of w = 10 -4 . The other two mocks were generated with 
a radial selection function using two exponential decaying models 
of completeness w (see Fig.Q} emulating apparent magnitude limit 
effects (see radial selection function in Kita ura et alJl2009h . The 
final mock galaxy samples have 350961 and 123679 particles us- 
ing the softer and steeper decaying selection functions respectively. 
The observer was set at the center of the box, i.e. at coordinates: 
X=250 Mpc/h, Y=250 Mpc/h, and Z=250 Mpc/h. 

The left panel on Fig. [2] shows a cell-to-cell comparison be- 
tween the dark matter distribution of the Millennium run (~ 10 10 
particles) gridded on a 256 3 mesh and a sub sample of 1 Million 
homo geneously selected mock galaxies of the lDe Lucia & Blaizotl 
catalogue. One can clearly see a deviation of the pixels with 
respect to the perfect slope of 45° . This effect is due to galaxy bias- 
ing. The right panel shows the analougous comparison with a Pois- 
son sampling using a homogeneous completeness of w = 10 -4 
which leaves about 1 Million particles. Here, we see a nearly per- 
fect scatter around the 45° slope demonstrating that our mocks do 
not include biasing. 



4.3 MAP results 

Here we calculate the maximum a posteriori solutions which we 
derived in the previous theoretical sections. There we assumed two 
different prior distributions for the matter field: a Gaussian and a 
lognormal prior. 



4.3. 1 Gaussian prior and Poissonian likelihood 

In this subsection we present the results given by the maximum 
a posteriori solution assuming a Gaussian prior and a Poissonian 
likelihood. The solution of Eqn. Q21 leads to a matter field which 
dramatically underestimates large overdensities (see Fig. O. This 
shows that the Gaussian prior cannot fit the underlying matter field 
which has a clearly non-Gaussian distribution with a minimum 
overdensity of ~ —1 up to maximal overdensities of about 1500 at 
the resolution we are looking at (~ 2 Mpc/h cell side length). The 
density peaks are highly suppressed with a Gaussian prior. This ef- 
fect is known from the Wiener filter as traditionally applied where 
the noise covarianc e is dependent on the signal (see discussion in 
iKitaura et alj2009h . Note, that the filter we are using here is a more 
accurate being based on the full Poissonian distribution and not 
only on the second order term as in the Wiener filter. 

4.3.2 Lognormal prior and Poissonian likelihood 

Here we present the results of the maximum a posteriori solution 
assuming a lognormal prior and a Poissonian likelihood. For that 
we solve the MAP Eqn. [32] 

We show in Fig.|4]the performance of the Poisson-lognormal 
filter with a homogeneous completeness. Panel a in Fig. [4] shows 
a slice through the matter distribution from the Millennium run. 
Panel b shows the mock galaxy sample. Panels c and d show the 
Poisson-lognormal filter and the LSQ filter reconstruction respec- 
tively. The performance depicted in cell-to-cell correlation plots 
shown in panels e and f demonstrate the superior behaviour of the 
Poisson-lognormal filter reconstruction in terms of higher corre- 
lation, smaller Euclidean distances and better alignment along the 
perfect correlation slope. The Poisson-lognormal filter recovers the 
density field up to overdensities above 1500 whereas the LSQ filter 
tends to underestimate the density field. 

We study the inhomogeneous completeness effects by select- 
ing dark matter particle subsamples with two different radial se- 
lection functions depicted in Fig. Q] In the upper panel of Fig. \5\ 
the inverse weighting scheme is shown to overestimate the density 
at low completeness (at the borders and corners of the cube) . This 
is in agreement with tests performed bv lKitaura etafl d2009h . The 
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Figure 4. Panel a: slice through the Millenium run dark matter particle simulation. Panel b: mock galaxy distribution with 10 6 particles. Panel c: reconstruction 
with the lognormal filter. Panel d: reconstruction with the LSQ-Wiener filter. Panel e: cell-to-cell correlation between the overdensity of the full three- 
dimensional reconstruction with the Lognormal filter (panel c) and the matter field (panel a). Panel f: cell-to-cell correlation between the full three-dimensional 
overdensity of the reconstruction with the LSQ-Wiener filter (panel d) and the matter field (panel a). The cell-to-cell correlation between the mock galaxy 
distribution (panel b) and the dark matter distribution can be seen on the right panel of Fig. [2 The plots were produced by calculating the mean over 15 
neighboring slices around slice 218 (Y~ 176 Mpc/h) through a 500 Mpc/h cube box with a 256 3 grid. 
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Figure 5. Panel a: inverse weigted mock galaxy distribution after applying a radial selection function (wi) leaving 350961 galaxies. Panels c: LSQ-Wiener 
filter reconstruction. Panels e: Lognormal filter reconstruction. The plots were produced by calculating the mean over 15 neighboring slices around slice 218 
(Y~ 176 Mpc/h) through a 500 Mpc/h cube box with a 256 3 grid. The performance depicted in cell-to-cell correlation plots are shown in the right hand side 
panels b, e and f. 
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Figure 7. Matter field reconstructions with the lognormal filter on a grid mesh with 128 3 and 64 3 cells for a uniform selection using about 10 6 mock galaxies. 
Panel a: mean over 9 slices through the reconstruction on a mesh with 128 3 cells around slice 109 (Y ~ 179 Mpc/h). Panel b: mean over 5 slices through 
the reconstruction on a mesh with 64 3 cells around slice 55 (Y ~ 179 Mpc/h). Panels c and d show the cell-to-cell statistics corresponding to the full 
reconstructions shown in panels a and b, respectively. 



LSQ filter, on the other hand, smooths the density more strongly 
in at low completeness regions and leads to a significantly lower 
Euclidean distance. The correlation coefficient is lower since the 
LSQ filtering suppresses the signal and gives a smooth version of 
the density field which is valid on larger scales (see panels c and d 
of Fig. |5]l, but does not reproduce small-scale features. The lower 
panels show the results coming from the Poisson-lognormal filter 
reconstruction. The density at low completeness is suppressed to 
zero due to the mean field used for this calculation (see Eq.[l6]). In 
regions of very low completeness the filter tends to favor the mean 
density. The statistical correlation shows to be clearly superior to 
the previous cases and the Euclidean distance with respect to the 
underlying matter field is far smaller (see panel f). The cell-to-cell 
correlation plot shows a scatter around the 45 slope and reproduces 
even the highest overdensities like the one at ^1600 which can also 
be seen in panel b. We perform the analogous study with the steeper 
radial selection function (see Fig.Q}. The results are shown in Fig.[6] 
and are consistent with the previously discussed ones. 

We perform the same study for two more resolutions: a mesh 
with 128 3 cells and a mesh with 64 3 cells for the same comov- 



ing box. First we grid the dark matter field coming from the Mil- 
lenium run on the lower resolution mesh and then we apply the 
radial selection w(r) using w — 10 -4 , w(r) — w\(r) and 
w(r) — W2(r). The results of the Poisson-lognormal filter recon- 
structions are shown in Figs. |7] and [8] We see a clear tendency to 
better recover the underlying matter field when using a lower reso- 
lution (compare Fig. |4] with Fig.|7]and Figs. [5] [6] with Fig. [8}. 

4.3.3 Matter statistics 

Finally, we calculate the matter statistics and the corresponding 
skewness and kurtosis for the dark matter field and the lognormal 
reconstructions corresponding to the incomplete mocks with selec- 
tion functions w± and wi (for particular expression s to calculate 
the m atter statistics, the skewness and the kurosis see lKitaura et al.l 
2009). The matter statistics represented in Fig. [9] shows consistent 
results for different grid resolutions (compare left, middle and right 
panels). After convolving the matter fields with a Gaussian kernel 
using a smoothing radius of 10 Mpc/h the matter distribution ap- 
pears to be closely lognormal distributed for all resolutions (see 
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Figure 8. Matter field reconstructions with the lognormal filter on a grid mesh with 128 3 and 64 3 cells for both w\ and W2 selection criteria. Panel a: same 
as panel a in previous figure for the case w±. Panel b: same as panel a for the case W2. Panels c and d show the cell-to-cell statistics corresponding to the full 
reconstructions shown in panels a and b, respectively. Panel e: same as panel a on a mesh with 64 3 . Panel f : same as panel e for the case of W2 . Panels g and h 
show the cell-to-cell statistics corresponding to the full reconstructions shown in panels e and f, respectively. 
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Figure 9. Matter statistics for the dark matter field from the Millenium run (black curve, red error bars) using about ~ 10 10 particles and the corresponding 
reconstructions using the selected mocks with radial completeness w\ (dashed curve, green error bars) and W2 (dashed-dotted curve, cyan error bars) having 
about ~ 10 5 particles for different resolutions (256 3 : left panels, 128 3 : middle panels and 64 3 : right panels). Upper panels (a, b and c): without smoothing. 
Lower panels: after convolution with a Gaussian kernel with smoothing radii of 5 Mpc/h (panels d, e and f) and 10 Mpc/h (panels g, h and i). The number of 
cells was counted for a logarithmic density binning of 0.2 in ln(l + 5m) for all cases except for panel a for which a binning of 0.4 was used. Also skewness 
(s, si and S2) and kurtosis (k, k\ and k2) in ln(l + 5m) are shown corresponding to the matter field, the reconstructions for the case w\ and the case W2, 
respectively. The error bars are given by the shot noise caused by the number counts of cells in each density bin without taking into account the uncertainties 
introduced by the completeness or the reconstruction method itself. 



panels at the bottom). The skewness (s, s± and S2) and the kurtosis 
(k, ki and fo) show some deviation from zero particularly in the 
reconstructed fields (skewness and kurtosis without subindex corre- 
spond to the dark matter field, with the subindex "1" to the selected 
sample with w± and with the subindex "2" to the selected sample 
with W2). However, their values are small which means that the 
distributions are not especially peaked or have significantly longer 
tails with respect to the lo gnormal distributio n. This result is con- 
sistent with observations ( Kita ura et al . 2009) where for a similar 
smoothing radius the matter field obtained from the Sloan Digital 



Sky survey (data release 6) was found to be close to lognormal dis- 
tributed. When convolving with a Gaussian kernel with a 5 Mpc/h 
smoothing radius (panels d, e and f) the distribution shows a tail 
towards larger densities with higher skewness and kurtosis than for 
the panels at the bottom which cannot be attributed to the uncer- 
tainty at the high densities shown by the large error bars caused by 
the low number counts in that regime. This deviation from the log- 
normal distribution is even better demonstrated in the upper panels 
which show the matter statistics without any additional smoothing. 
The results show that the multivariate lognormal prior distribution 
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does not impose a lognormal matter field statistics to the recovered 
density field. Here the prior is subdominant with respect to the data 
similarly to the case of the LSQ-Wiener reconstruction which can 
lead to non-Ga ussian statistics ev en though it is based on a Gaus- 
sian prior (see lKitaura et 

However, we also observe several effects causing a deviation 
in the reconstructed fields from the true matter field statistics. The 
reconstructions tend to overestimate the number of cells around the 
mean density (see the peaks of the distribution in the upper pan- 
els of Fig. [9). This trend is more acute for the stronger sampled 
mock for which W2 was used. This can be seen by comparing the 
dashed-dotted curves with the dashed curves and the kurtosis k\ 
with (fe is always larger than hi). We also observe that the re- 
constructions underestimate the number of cells in the extremely 
underdense regions (8 < — 0.6). In addition, we investigated the 
statistics of the reconstructed fields based on the homogeneously 
sampled mock data using w = 10 -4 and found that the densities 
are distributed in a very similar way to the reconstructed matter 
fields with w\. These effects are caused by the conservative char- 
acter of the reconstruction method. The maximum a posteriori solu- 
tion leads to a stronger smoothing in the undersampled low-density 
regions and produces a larger number of cells with densities closer 
to the mean. 



5 CONCLUSIONS 

In this work we have presented a general expression for the 
Poisson-lognormal filter given an arbitrary nonlinear galaxy bias. 
We derived this filter as the maximum a posteriori solution assum- 
ing a lognormal prior distribution for the matter field with a con- 
stant mean field and modeling the observed galaxy distribution by 
a Poissonian process (see Eq.|26l). 

We have performed a three-dimensional implementation of 
this filter with a very efficient Newton-Krylov inversion scheme 
(see section [3]). Furthermore, we have tested it for a linear galaxy 
bias relation and compared the results with other density field esti- 
mators commonly used in the literature (e.g. the inverse weighting 
scheme and the least squares (LSQ) Wiener filter (see section |4])). 

We also found that the solution of Eqn. Q21 assuming a Gaus- 
sian prior distribution for the matter field, leads to a reconstruc- 
tion which clearly underestimates large overdensities (see Fig. [3]). 
This shows that the Gaussian prior cannot fit the underlying matter 
field which has a clearly non-Gaussian distribution with a mini- 
mum overdensity of 5 ~ — 1 up to maximal overdensities of about 
S ~ 1700 for a resolution of ~ 2 Mpc/h. The density peaks are 
highly suppressed with the Gaussian prior. This effect is known 
from the Wiener filter as traditionally applied in whic h the noise co- 
varian ce is dependent on the signal (see discussion in lKitaura et al.1 
l2009h . 

However, we have seen that even the LSQ-Wiener filter fails 
for high overdensities (5m > 100). We showed in appendix [A] 
that the LSQ-filter is the optimal linear filter under a Poisson 
noise assumption and does not neglect any signal to noise corre- 
lation, c ontrary to what has been assumed in literature (see for 
example Izaroubi et all II 19951 : ISeliakl 1 19981 : lErdogdu et all 120041 : 
iKitaura & EnBlinll2008h . The LSQ filter is the optimal linear fil- 
ter only up to second order statistics and thus is less well suited 
to distributions with high skewness or long tails such as the log- 
normal distribution. Another reason for the inferior performance of 
the LSQ-filter with respect to the Poisson-lognormal filter is its 



linearity. Note, that the relation at overdensities <5 ^> 1 is highly 
nonlinear. 

The one-dimensional lognormal probability distribution is 
known to fit well the matter distribution up to overdensities of about 
5 ^100 as found bv lKavo etall fcoOD. Our results show, however, 
good agreement for overdensities even above 5 ^1000 which ex- 
ceeds by one order of magnitude the regime in which the lognor- 
mal is expected to be valid. This is because in our filter the log- 
normal assumption enters as a prior distribution function, but the 
maximum a posteriori solution i s also conditioned on the data. In a 
similar way Kit aura et ai] j2009h was able to recover a highly non- 
Gaussian distributed matter field from the SDSS dr6 after using 
LSQ-Wiener filtering which according to the Bayesian formalism 
assumes a Gaussian distribution. Assuming that the galaxy bias is 
known we find that the Poisson-lognormal filter is able to recover 
the matter density fields down to scales of about > 2 Mpc/h. How- 
ever, our study of the matter statistics comparing the dark matter 
with the reconstructed fields shows that the Poisson-lognormal fil- 
ter fails to recover underdense regions for 8 < — 0.6. At lower den- 
sities the recovered field is smoothed out due to the conservative 
maximum a posteriori solution. 

Our work shows a great improvement with respect to previous 
filters in recovering the matter density field from a point source dis- 
tribution. Still much work has to be done to further analyse the sta- 
tistical properties of the cosmological structure. Nevertheless, the 
nonlinear reconstruction method derived in this work could be of 
great interest for large scale structure density field reconstructions 
taking a galaxy distribution or even some other observables like the 
Lyman alpha forest. 
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APPENDIX A: THE LSQ FILTER 

Here we show a derivation of the LSQ filter which does not require 
a data degradation model with an additive noise term. Let us adopt 
here the usual notation for the data: d = 5°. The data vector is 
accordingly defined by: 



N° 



Wl, 



(Al) 



for a definition of the mask operator W see section lT. 1.11 We define 
here the signal vector s as the matter overdensity field: s = 5m. In 
the linear approximation we try to find a filter F which applied to 
the data d gives an estimate of the signal s of the form: 



(s)lsq Fd. 



(A2) 



This filter should mini mize the follow i ng quantity in the least 
squares approach (see IWienerl Il949l : [Rybicki & P ressl Il992l : 
IZaroubi et all 19951) : 

A = ((Fd-s) 2 ) 

= F(dd t )F t - F(ds^) - (sd^F 1 " + (ss f ). (A3) 

As lKitaura & EnBlinl (l2008h pointed out it is important to note 
that the ensemble average ({ }) goes over the galaxy and matter 
field realizations and the filter is thus different from the Wiener 
filter as derived in a Bayesian framework. We define here the 
global ensemble average by: ({ }) = (({}) g )M- Here ({ }) g = 
({ })(iV||A o ) = X}iv°= Ppois(N° | w\){ } denotes an ensemble 
average over the Poissonian distribution with the expected number 
of galaxy counts given by the Poissonian ensemble average: A° = 
w\= (iV g ) g ,and({})M = ({})(s M \p M ) = f d5 M P(S M | p M ) 
being the ensemble average over all possible matter density real- 
izations with some prior distribution P(Sm | P m ) w * m Pm being a 
set of parameters which determine the matter field, say the cosmo- 
logical parameters. We impose (5m) m — 0. 

Recalling t he d eri vations done by IWiened dl949h : 
iRvbicki & Pressl dl992h : IZaroubi et all dl995h we find mini- 
mizing the action with respect to the filter: 



OA 
OF 



: 0, 



the following LSQ filter expression: 

F= (ad t )(dd t )- 1 . 



(A4) 



(A5) 



Traditionally one would then define a data degradation model with 
an additive noise term of the form: d — Rs + e, with R being some 
response operator. Then substituting this data model in Eq. lA5l and 
neglecting noise to signal correlat ion terms one wou ld obtain a final 
expression for the LSQ filter (see Zaroubi eTaD 19951) . 



Al Signal to noise correlation 

One can show that the noise is actually uncorrected with the signal 
by making the following definition: 

e° = iV° t - (iV g ° ;i ) g , (A6) 

and then calculating the correlation: 

fe°(Ki)«>« = (KA N li)z - (Nli)&(K,j) g }& = o- (A7) 
Note, that this implies: (e°(£g,j)g) g = and thus also (es^) g = 0. 
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A2 LSQ filter derivation without the additive noise 
assumption 

However, one does not even need to use the additive noise assump- 
tion to derive the LSQ filter. Let us show here how to make such a 
derivation. We define the observed galaxy number counts per cell i 
as: 

Nl t = N g (w t + 6l z ). (A8) 

The corresponding ensemble average over all possible galaxy real- 
izations is: 

(NZJ s = N g Wi(l + 6 g ,i). (A9) 
Recalling the linear bias relation: 

S s ,i = ^2b id S h (A10) 

3 

we can then calculate with the above definitions the signal to data 
correlation matrix: 

(sidj) = (S M ,iSg,j) = {{SM,iSg,j) s )M 

= (<$M,z(£g,j)g)M = wj y^bjjf (Sm^Smj^m- (All) 
j' 

We also have to calculate the data autocorrelation matrix: 



(didj) = (6g,i&&,j) = (( 



N, 



Nt 



-WiWj. (A12) 



Here we need a model for the two-point number count statistics. 
Note, that we can introduce here Poissonity: 

(N^N^), = (Nl t ) s (Nl J ) s + (N^)^. (A13) 

With the additional matter field ensemble average we get: 

((iVg,i>g(A^g ;J -> g >M = N 2 g WiWj I 1 + y^ j b iik y^ j bj i i(5M,kSM,i)< 

\ k i 



We can define the noise covariance matrix as: 

Ni,j ee i<(iV g ,^> g -(iV g , i > g (iV g J > g ) r 



h3 — 2 

The LSQ filter can be then written as: 



(A14) 



(A15) 



(A16) 



Ne 3 ,J 



X I W 3' b 3',k ^ b j,k> (5 M ,kSM,k>)MWj + 
\ k k' 

The corresponding matrix notation of the LSQ filter yields: 

F = SR 1 " (RSR 1 " + n) _1 , (A17) 

with S = (<5m<5] vI )m and R = WB (see section fl~ 1.1 1 for a defi- 
nition of the bias operator B). This data-space expression is equiv- 
alent to the s ignal- space representati on (for a demonstration see 
appendix C in lKitaura & EnBlinll2008h : 



F = (V 1 + R+N" 1 ^ 1 R f N _1 . 



We conclude that the LSQ filter is the optimal linear filter under a 
Poisson noise assumption. We have shown that this filter does not 
neglect any signal to noise correlation. 



APPENDIX B: LOGNORMAL PRIOR AND GAUSSIAN 
LIKELIHOOD 

For completeness we derive a nonlinear filter which assumes alog- 
normal prior and a Gaussian likelihood. Following Sheth (1995) 
one could use the Wiener filter with a data transformation and ap- 
ply it to recover non-Gaussian distributed fields. The problem in 
such a model is that one requires a multiplicative noise assumption 
of the form: 

d'i = 5i€i + e'i 
ln(^) = ln(l + fc)+ln(ej) 

d% — Si -\~ €i . 



(Bl) 



for each cell i, with d% = ln(c^), Si = ln(l + Si) and a = ln(e^). 
Note, that with such a data model one could easily apply the Wiener 
filter assuming that the signal s and the noise e are Gaussian dis- 
tributed. 



Bl Additive noise model 

However, one may rather prefer a data model with an additive noise 
term as commonly us ed in the literature (see e.g. IZaroubi et al.l 
1995; Tegmark 1997). We define therefore a data model of the 
form: 



d = RS + e, 
including in the signal higher order terms: 
S — exp(s) — 1. 



(B2) 



(B3) 



One can then assume the signal to be lognormal distributed and the 
noise to be Gaussian distributed and signal-independent. 



B2 Gaussian likelihood 

Let us write the log-likelihood as: 

ln£oc -~ ^e t N" 1 e-ln(det(N))) + c. (B4) 
Making the substitution e = d — RS we get: 

e+N^e = d^N^d + ^R+N^Rd - 5 t R t N" 1 - d^N^RS. 

(B5) 

To find the maximum a posteriori solution we have to calculate the 
derivative of the likelihood with respect to the signal s: 



Ed In d _ ^ ® ^ n 



ds k 



From Eq. lB3l we get: 



ds k 



■ exp(si)8f; k . 



(B6) 



(B7) 



Assuming a signal independent noise yields: 
d\n£i 



(A18) 



jlm 
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Poisson-Lognormal Reconstruction 

Combining these results with the derivative of the lognormal prior 
(Eq. [19]) leads to: 

E^fe-/^)= (B9) 

3 

\jlm jl I 
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